﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices;

namespace Rmath64Dll
{
    public class Rmath
    {
        private const double UNDEFINE = -99999;
        /* Names from Solaris math.h */
        public const double M_E = 2.718281828459045235360287471353;	/* e */
        public const double M_LOG2E = 1.442695040888963407359924681002;	/* log2(e) */
        public const double M_LOG10E = 0.434294481903251827651128918917;	/* log10(e) */
        public const double M_LN2 = 0.693147180559945309417232121458;	/* ln(2) */
        public const double M_LN10 = 2.302585092994045684017991454684;/* ln(10) */
        public const double M_PI = 3.141592653589793238462643383280;	/* pi */
        public const double M_2PI = 6.283185307179586476925286766559;	/* 2*pi */
        public const double M_PI_2 = 1.570796326794896619231321691640;	/* pi/2 */
        public const double M_PI_4 = 0.785398163397448309615660845820;	/* pi/4 */
        public const double M_1_PI = 0.318309886183790671537767526745;/* 1/pi */
        public const double M_2_PI = 0.636619772367581343075535053490;	/* 2/pi */
        public const double M_2_SQRTPI = 1.128379167095512573896158903122;/* 2/sqrt(pi) */
        public const double M_SQRT2 = 1.414213562373095048801688724210;	/* sqrt(2) */
        public const double M_SQRT1_2 = 0.707106781186547524400844362105;	/* 1/sqrt(2) */

        /* R-Specific Constants */
        public const double M_SQRT_3 = 1.732050807568877293527446341506;/* sqrt(3) */
        public const double M_SQRT_32 = 5.656854249492380195206754896838;	/* sqrt(32) */
        public const double M_LOG10_2 = 0.301029995663981195213738894724;	/* log10(2) */
        public const double M_SQRT_PI = 1.772453850905516027298167483341;/* sqrt(pi) */
        public const double M_1_SQRT_2PI = 0.398942280401432677939946059934;	/* 1/sqrt(2pi) */
        public const double M_SQRT_2dPI = 0.797884560802865355879892119869;/* sqrt(2/pi) */
        public const double M_LN_SQRT_PI = 0.572364942924700087071713675677;	/* log(sqrt(pi))== log(pi)/2 */
        public const double M_LN_SQRT_2PI = 0.918938533204672741780329736406;	/* log(sqrt(2*pi)) == log(2*pi)/2 */
        public const double M_LN_SQRT_PId2 = 0.225791352644727432363097614947;/* log(sqrt(pi/2)) */

        #region Rmath.dll functions----------------------------------------------------------------------------

        [DllImport(@"Rmath64.dll", EntryPoint = "R_pow")]
        private static extern double R_pow_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "R_pow_di")]
        private static extern double R_pow_di_dll(double a, int b);

        [DllImport(@"Rmath64.dll", EntryPoint = "norm_rand")]
        private static extern double norm_rand_dll();

        [DllImport(@"Rmath64.dll", EntryPoint = "unif_rand")]
        private static extern double unif_rand_dll();

        [DllImport(@"Rmath64.dll", EntryPoint = "exp_rand")]
        private static extern double exp_rand_dll();

        [DllImport(@"Rmath64.dll", EntryPoint = "set_seed")]
        private static extern void set_seed_dll(uint a, uint b);

        [DllImport(@"Rmath64.dll", EntryPoint = "get_seed")]
        private static extern void get_seed_dll(uint[] a, uint[] b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnorm4")]
        private static extern double dnorm_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnorm5")]
        private static extern double pnorm_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnorm5")]
        private static extern double qnorm_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rnorm")]
        private static extern double rnorm_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnorm_both")]
        private static extern void pnorm_both_dll(double a, double[] c, double[] d, int e, int f);/* both tails */

        [DllImport(@"Rmath64.dll", EntryPoint = "dunif")]
        private static extern double dunif_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "punif")]
        private static extern double punif_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qunif")]
        private static extern double qunif_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "runif")]
        private static extern double runif_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dgamma")]
        private static extern double dgamma_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pgamma")]
        private static extern double pgamma_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qgamma")]
        private static extern double qgamma_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rgamma")]
        private static extern double rgamma_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dbeta")]
        private static extern double dbeta_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pbeta")]
        private static extern double pbeta_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qbeta")]
        private static extern double qbeta_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rbeta")]
        private static extern double rbeta_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dlnorm")]
        private static extern double dlnorm_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "plnorm")]
        private static extern double plnorm_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qlnorm")]
        private static extern double qlnorm_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rlnorm")]
        private static extern double rlnorm_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dchisq")]
        private static extern double dchisq_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "pchisq")]
        private static extern double pchisq_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qchisq")]
        private static extern double qchisq_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rchisq")]
        private static extern double rchisq_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnchisq")]
        private static extern double dnchisq_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnchisq")]
        private static extern double pnchisq_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnchisq")]
        private static extern double qnchisq_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rnchisq")]
        private static extern double rnchisq_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "df")]
        private static extern double df_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pf")]
        private static extern double pf_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qf")]
        private static extern double qf_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rf")]
        private static extern double rf_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dt")]
        private static extern double dt_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "pt")]
        private static extern double pt_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qt")]
        private static extern double qt_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rt")]
        private static extern double rt_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "dbinom")]
        private static extern double dbinom_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pbinom")]
        private static extern double pbinom_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qbinom")]
        private static extern double qbinom_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rbinom")]
        private static extern double rbinom_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "rmultinom")]
        private static extern void rmultinom_dll(int a, double[] b, int c, int[] d);

        [DllImport(@"Rmath64.dll", EntryPoint = "dcauchy")]
        private static extern double dcauchy_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pcauchy")]
        private static extern double pcauchy_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qcauchy")]
        private static extern double qcauchy_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rcauchy")]
        private static extern double rcauchy_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dexp")]
        private static extern double dexp_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "pexp")]
        private static extern double pexp_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qexp")]
        private static extern double qexp_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rexp")]
        private static extern double rexp_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "dgeom")]
        private static extern double dgeom_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "pgeom")]
        private static extern double pgeom_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qgeom")]
        private static extern double qgeom_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rgeom")]
        private static extern double rgeom_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "dhyper")]
        private static extern double dhyper_dll(double a, double b, double c, double d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "phyper")]
        private static extern double phyper_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "qhyper")]
        private static extern double qhyper_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "rhyper")]
        private static extern double rhyper_dll(double a, double b, double c);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnbinom")]
        private static extern double dnbinom_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnbinom")]
        private static extern double pnbinom_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnbinom")]
        private static extern double qnbinom_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rnbinom")]
        private static extern double rnbinom_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnbinom_mu")]
        private static extern double dnbinom_mu_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnbinom_mu")]
        private static extern double pnbinom_mu_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnbinom_mu")]
        private static extern double qnbinom_mu_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rnbinom_mu")]
        private static extern double rnbinom_mu_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dpois")]
        private static extern double dpois_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "ppois")]
        private static extern double ppois_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qpois")]
        private static extern double qpois_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rpois")]
        private static extern double rpois_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "dweibull")]
        private static extern double dweibull_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pweibull")]
        private static extern double pweibull_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qweibull")]
        private static extern double qweibull_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rweibull")]
        private static extern double rweibull_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dlogis")]
        private static extern double dlogis_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "plogis")]
        private static extern double plogis_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qlogis")]
        private static extern double qlogis_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "rlogis")]
        private static extern double rlogis_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnbeta")]
        private static extern double dnbeta_dll(double a, double b, double c, double d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnbeta")]
        private static extern double pnbeta_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnbeta")]
        private static extern double qnbeta_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "rnbeta")]
        private static extern double rnbeta_dll(double a, double b, double c);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnf")]
        private static extern double dnf_dll(double a, double b, double c, double d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnf")]
        private static extern double pnf_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnf")]
        private static extern double qnf_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "dnt")]
        private static extern double dnt_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pnt")]
        private static extern double pnt_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "qnt")]
        private static extern double qnt_dll(double a, double b, double c, int d, int e);

        [DllImport(@"Rmath64.dll", EntryPoint = "ptukey")]
        private static extern double ptukey_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "qtukey")]
        private static extern double qtukey_dll(double a, double b, double c, double d, int e, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "dwilcox")]
        private static extern double dwilcox_dll(double a, double b, double c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "pwilcox")]
        private static extern double pwilcox_dll(double a, double b, double c, int d, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "qwilcox")]
        private static extern double qwilcox_dll(double a, double b, double c, int d, int f);

        [DllImport(@"Rmath64.dll", EntryPoint = "rwilcox")]
        private static extern double rwilcox_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dsignrank")]
        private static extern double dsignrank_dll(double a, double b, int c);

        [DllImport(@"Rmath64.dll", EntryPoint = "psignrank")]
        private static extern double psignrank_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "qsignrank")]
        private static extern double qsignrank_dll(double a, double b, int c, int d);

        [DllImport(@"Rmath64.dll", EntryPoint = "rsignrank")]
        private static extern double rsignrank_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "gammafn")]
        private static extern double gammafn_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "lgammafn")]
        private static extern double lgammafn_dll(double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "lgammafn_sign")]
        private static extern double lgammafn_sign_dll(double a, int[] b);

        [DllImport(@"Rmath64.dll", EntryPoint = "dpsifn")]
        private static extern void dpsifn_dll(double a, int b, int c, int d, double[] e, int[] f, int[] g);

        [DllImport(@"Rmath64.dll", EntryPoint = "psigamma")]
        private static extern double psigamma_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "digamma")]
        private static extern double digamma_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "trigamma")]
        private static extern double trigamma_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "tetragamma")]
        private static extern double tetragamma_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "pentagamma")]
        private static extern double pentagamma_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "beta")]
        private static extern double beta_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "lbeta")]
        private static extern double lbeta_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "choose")]
        private static extern double choose_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "lchoose")]
        private static extern double lchoose_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_i")]
        private static extern double bessel_i_dll(double a, double b, double c);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_j")]
        private static extern double bessel_j_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_k")]
        private static extern double bessel_k_dll(double a, double b, double c);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_y")]
        private static extern double bessel_y_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_i_ex")]
        private static extern double bessel_i_ex_dll(double a, double b, double c, double[] d);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_j_ex")]
        private static extern double bessel_j_ex_dll(double a, double b, double[] c);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_k_ex")]
        private static extern double bessel_k_ex_dll(double a, double b, double c, double[] d);

        [DllImport(@"Rmath64.dll", EntryPoint = "bessel_y_ex")]
        private static extern double bessel_y_ex_dll(double a, double b, double[] c);

        [DllImport(@"Rmath64.dll", EntryPoint = "hypot")]
        private static extern double hypot_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "pythag")]
        private static extern double pythag_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "expm1")]
        private static extern double expm1_dll(double a); /* = exp(x)-1 {care for small x} */

        [DllImport(@"Rmath64.dll", EntryPoint = "log1p")]
        private static extern double log1p_dll(double a); /* = log(1+x) {care for small x} */

        [DllImport(@"Rmath64.dll", EntryPoint = "imax2")]
        private static extern int imax2_dll(int a, int b);

        [DllImport(@"Rmath64.dll", EntryPoint = "imin2")]
        private static extern int imin2_dll(int a, int b);

        [DllImport(@"Rmath64.dll", EntryPoint = "fmax2")]
        private static extern double fmax2_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "fmin2")]
        private static extern double fmin2_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "sign")]
        private static extern double sign_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "fprec")]
        private static extern double fprec_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "fround")]
        private static extern double fround_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "fsign")]
        private static extern double fsign_dll(double a, double b);

        [DllImport(@"Rmath64.dll", EntryPoint = "ftrunc")]
        private static extern double ftrunc_dll(double a);

        [DllImport(@"Rmath64.dll", EntryPoint = "log1pmx")]
        private static extern double log1pmx_dll(double a); /* Accurate log(1+x) - x, {care for small x} */

        [DllImport(@"Rmath64.dll", EntryPoint = "lgamma1p")]
        private static extern double lgamma1p_dll(double a);/* accurate log(gamma(x+1)), small x (0 < x < 0.5) */

        [DllImport(@"Rmath64.dll", EntryPoint = "logspace_add")]
        private static extern double logspace_add_dll(double logx, double logy);

        [DllImport(@"Rmath64.dll", EntryPoint = "logspace_sub")]
        private static extern double logspace_sub_dll(double logx, double logy);

        #endregion Rmath.dll functions----------------------------------------------------------------------------

        #region 添加wrapper------------------------------------------------------------------------------------
        public static void setseed(uint a, uint b)
        {
            set_seed_dll(a, b);
        }

        #region Normal Distribution----------------------------------------------------------------------
        public static double dnorm(double x, double mean = 0, double sd = 1, Boolean log = false)
        {
            return dnorm_dll(x, mean, sd, log ? 1 : 0);
        }
        public static double[] dnorm(double[] x, double[] mean = null, double[] sd = null, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnorm_dll(x[i], mean == null ? 0 : mean[i], sd == null ? 1 : sd[i], log ? 1 : 0);

            return result;
        }

        public static double pnorm(double x, double mean = 0, double sd = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return pnorm_dll(x, mean, sd, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pnorm(double[] x, double[] mean = null, double[] sd = null, Boolean lower_tail = true, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnorm_dll(x[i], mean == null ? 0 : mean[i], sd == null ? 1 : sd[i], lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qnorm(double x, double mean = 0, double sd = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qnorm_dll(x, mean, sd, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qnorm(double[] x, double[] mean = null, double[] sd = null, Boolean lower_tail = true, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnorm_dll(x[i], mean == null ? 0 : mean[i], sd == null ? 1 : sd[i], lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rnorm(int n, double mean = 0, double sd = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rnorm_dll(mean, sd);

            return result;
        }
        #endregion  Normal Distribution-----------------------------------------------------------------

        #region Uniform Distribution--------------------------------------------------------------
        public static double dunif(double x, double min = 0, double max = 1, Boolean log = false)
        {
            return dunif_dll(x, min, max, log ? 1 : 0);
        }
        public static double[] dunif(double[] x, double min = 0, double max = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dunif_dll(x[i], min, max, log ? 1 : 0);

            return result;
        }

        public static double punif(double q, double min = 0, double max = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return punif_dll(q, min, max, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] punif(double[] q, double min = 0, double max = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = punif_dll(q[i], min, max, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qunif(double p, double min = 0, double max = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qunif_dll(p, min, max, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qunif(double[] p, double min = 0, double max = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qunif_dll(p[i], min, max, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] runif(int n, double min = 0, double max = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = runif_dll(min, max);

            return result;
        }

        #endregion Uniform Distribution-------------------------------------------------------------

        #region Gamma Distribution---------------------------------------------------------
        public static double dgamma(double x, double shape, double rate = 1, double scale = UNDEFINE, Boolean log = false)
        {
            return dgamma_dll(x, shape, scale == UNDEFINE ? 1 / rate : scale, log ? 1 : 0);
        }
        public static double[] dgamma(double[] x, double shape, double rate = 1, double scale = UNDEFINE, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dgamma_dll(x[i], shape, scale == UNDEFINE ? 1 / rate : scale, log ? 1 : 0);

            return result;
        }

        public static double pgamma(double x, double shape, double rate = 1, double scale = UNDEFINE, Boolean lower_tail = true, Boolean log = false)
        {
            return pgamma_dll(x, shape, scale == UNDEFINE ? 1 / rate : scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pgamma(double[] x, double shape, double rate = 1, double scale = UNDEFINE, Boolean lower_tail = true, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pgamma_dll(x[i], shape, scale == UNDEFINE ? 1 / rate : scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qgamma(double x, double shape, double rate = 1, double scale = UNDEFINE, Boolean lower_tail = true, Boolean log = false)
        {
            return qgamma_dll(x, shape, scale == UNDEFINE ? 1 / rate : scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qgamma(double[] x, double shape, double rate = 1, double scale = UNDEFINE, Boolean lower_tail = true, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qgamma_dll(x[i], shape, scale == UNDEFINE ? 1 / rate : scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rgamma(int n, double shape, double rate = 1, double scale = UNDEFINE)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rgamma_dll(shape, scale == UNDEFINE ? 1 / rate : scale);

            return result;
        }
        #endregion Gamma Distribution------------------------------------------------------

        #region Beta Distribution  and Non-central Beta Distribution--------------------------
        public static double dbeta(double x, double shape1, double shape2, Boolean log = false)
        {
            return dbeta_dll(x, shape1, shape2, log ? 1 : 0);
        }
        public static double dbeta(double x, double shape1, double shape2, double ncp, Boolean log = false)
        {
            return dnbeta_dll(x, shape1, shape2, ncp, log ? 1 : 0);
        }
        public static double[] dbeta(double[] x, double shape1, double shape2, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dbeta_dll(x[i], shape1, shape2, log ? 1 : 0);

            return result;
        }
        public static double[] dbeta(double[] x, double shape1, double shape2, double ncp, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnbeta_dll(x[i], shape1, shape2, ncp, log ? 1 : 0);

            return result;
        }

        public static double pbeta(double q, double shape1, double shape2, Boolean lower_tail = true, Boolean log = false)
        {
            return pbeta_dll(q, shape1, shape2, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double pbeta(double q, double shape1, double shape2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return pnbeta_dll(q, shape1, shape2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pbeta(double[] q, double shape1, double shape2, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pbeta_dll(q[i], shape1, shape2, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] pbeta(double[] q, double shape1, double shape2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnbeta_dll(q[i], shape1, shape2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qbeta(double p, double shape1, double shape2, Boolean lower_tail = true, Boolean log = false)
        {
            return qbeta_dll(p, shape1, shape2, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double qbeta(double p, double shape1, double shape2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return qnbeta_dll(p, shape1, shape2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qbeta(double[] p, double shape1, double shape2, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qbeta_dll(p[i], shape1, shape2, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] qbeta(double[] p, double shape1, double shape2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnbeta_dll(p[i], shape1, shape2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rbeta(int n, double shape1, double shape2, double ncp = 0)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            if (ncp == 0)
                for (var i = 0; i < n; i++)
                    result[i] = rbeta_dll(shape1, shape2);
            else
            {
                var temp1 = rchisq(n, 2 * shape1, ncp);
                var temp2 = rchisq(n, 2 * shape2);
                for (var i = 0; i < n; i++)
                {
                    result[i] = temp1[i] / (temp1[i] + temp2[i]);
                }
            }
            return result;
        }
        #endregion Beta Distribution and Non-central Beta Distribution-------------------------

        #region Lognormal Distribution------------------------------------------------
        public static double dlnorm(double x, double meanlog = 0, double sdlog = 1, Boolean log = false)
        {
            return dlnorm_dll(x, meanlog, sdlog, log ? 1 : 0);
        }
        public static double[] dlnorm(double[] x, double meanlog = 0, double sdlog = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dlnorm_dll(x[i], meanlog, sdlog, log ? 1 : 0);

            return result;
        }

        public static double plnorm(double q, double meanlog = 0, double sdlog = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return plnorm_dll(q, meanlog, sdlog, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] plnorm(double[] q, double meanlog = 0, double sdlog = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = plnorm_dll(q[i], meanlog, sdlog, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qlnorm(double p, double meanlog = 0, double sdlog = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qlnorm_dll(p, meanlog, sdlog, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qlnorm(double[] p, double meanlog = 0, double sdlog = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qlnorm_dll(p[i], meanlog, sdlog, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rlnorm(int n, double meanlog = 0, double sdlog = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rlnorm_dll(meanlog, sdlog);
            return result;
        }
        #endregion Lognormal Distribution--------------------------------------------

        #region Chi-squared Distribution and Non-central Chi-squared Distribution----------------------------------------
        public static double dchisq(double x, double df, Boolean log = false)
        {
            return dchisq_dll(x, df, log ? 1 : 0);
        }
        public static double dchisq(double x, double df, double ncp, Boolean log = false)
        {
            return dnchisq_dll(x, df, ncp, log ? 1 : 0);
        }
        public static double[] dchisq(double[] x, double df, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dchisq_dll(x[i], df, log ? 1 : 0);

            return result;
        }
        public static double[] dchisq(double[] x, double df, double ncp, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnchisq_dll(x[i], df, ncp, log ? 1 : 0);

            return result;
        }

        public static double pchisq(double q, double df, Boolean lower_tail = true, Boolean log = false)
        {
            return pchisq_dll(q, df, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double pchisq(double q, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return pnchisq_dll(q, df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pchisq(double[] q, double df, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pchisq_dll(q[i], df, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] pchisq(double[] q, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnchisq_dll(q[i], df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qchisq(double p, double df, Boolean lower_tail = true, Boolean log = false)
        {
            return qchisq_dll(p, df, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double qchisq(double p, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return qnchisq_dll(p, df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qchisq(double[] p, double df, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qchisq_dll(p[i], df, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] qchisq(double[] p, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnchisq_dll(p[i], df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rchisq(int n, double df, double ncp = 0)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            if (ncp == 0)
                for (var i = 0; i < n; i++)
                    result[i] = rchisq_dll(df);
            else
            {
                for (var i = 0; i < n; i++)
                {
                    result[i] = rnchisq_dll(df, ncp);
                }
            }
            return result;
        }
        #endregion Chi-squared Distribution and Non-central Chi-squared Distribution-------------------------------------

        #region F Distibution and Non-central F Distribution-----------------------------------------------
        public static double df(double x, double df1, double df2, Boolean log = false)
        {
            return df_dll(x, df1, df2, log ? 1 : 0);
        }
        public static double df(double x, double df1, double df2, double ncp, Boolean log = false)
        {
            return dnf_dll(x, df1, df2, ncp, log ? 1 : 0);
        }
        public static double[] df(double[] x, double df1, double df2, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = df_dll(x[i], df1, df2, log ? 1 : 0);

            return result;
        }
        public static double[] df(double[] x, double df1, double df2, double ncp, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnf_dll(x[i], df1, df2, ncp, log ? 1 : 0);

            return result;
        }

        public static double pf(double q, double df1, double df2, Boolean lower_tail = true, Boolean log = false)
        {
            return pf_dll(q, df1, df2, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double pf(double q, double df1, double df2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return pnf_dll(q, df1, df2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pf(double[] q, double df1, double df2, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pf_dll(q[i], df1, df2, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] pf(double[] q, double df1, double df2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnf_dll(q[i], df1, df2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qf(double p, double df1, double df2, Boolean lower_tail = true, Boolean log = false)
        {
            return qf_dll(p, df1, df2, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double qf(double p, double df1, double df2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return qnf_dll(p, df1, df2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qf(double[] p, double df1, double df2, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qf_dll(p[i], df1, df2, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] qf(double[] p, double df1, double df2, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnf_dll(p[i], df1, df2, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rf(int n, double df1, double df2, double ncp = 0)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            if (ncp == 0)
                for (var i = 0; i < n; i++)
                    result[i] = rf_dll(df1, df2);
            else
            {
                var temp1 = rchisq(n, df1, ncp);
                var temp2 = rchisq(n, df2);
                for (var i = 0; i < n; i++)
                {
                    result[i] = (temp1[i] / df1) / (temp2[i] / df2);
                }
            }
            return result;
        }
        #endregion F Distibution and Non-central F Distribution-----------------------------------------------

        #region F Distibution and Non-central Student t Distribution-------------------------------
        public static double dt(double x, double df, Boolean log = false)
        {
            return dt_dll(x, df, log ? 1 : 0);
        }
        public static double dt(double x, double df, double ncp, Boolean log = false)
        {
            return dnt_dll(x, df, ncp, log ? 1 : 0);
        }
        public static double[] dt(double[] x, double df, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dt_dll(x[i], df, log ? 1 : 0);

            return result;
        }
        public static double[] dt(double[] x, double df, double ncp, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnt_dll(x[i], df, ncp, log ? 1 : 0);

            return result;
        }

        public static double pt(double q, double df, Boolean lower_tail = true, Boolean log = false)
        {
            return pt_dll(q, df, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double pt(double q, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return pnt_dll(q, df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pt(double[] q, double df, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pt_dll(q[i], df, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] pt(double[] q, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnt_dll(q[i], df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qt(double p, double df, Boolean lower_tail = true, Boolean log = false)
        {
            return qt_dll(p, df, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double qt(double p, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            return qnt_dll(p, df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qt(double[] p, double df, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qt_dll(p[i], df, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] qt(double[] p, double df, double ncp, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnt_dll(p[i], df, ncp, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rt(int n, double df)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rt_dll(df);
            return result;
        }
        public static double[] rt(int n, double df, double ncp = 0)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            var temp1 = rnorm(n, ncp);
            var temp2 = rchisq(n, df);
            for (var i = 0; i < n; i++)
            {
                result[i] = temp1[i] / Math.Sqrt(temp2[i] / df);
            }

            return result;
        }

        #endregion F Distibution and Non-central Student t Distribution-------------------------------

        #region Binomial Distribution-------------------------------------------------------
        public static double dbinom(double x, double size, double prob, Boolean log = false)
        {
            return dbinom_dll(x, size, prob, log ? 1 : 0);
        }
        public static double[] dbinom(double[] x, double size, double prob, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dbinom_dll(x[i], size, prob, log ? 1 : 0);

            return result;
        }

        public static double pbinom(double q, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return pbinom_dll(q, size, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pbinom(double[] q, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pbinom_dll(q[i], size, prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qbinom(double p, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return qbinom_dll(p, size, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qbinom(double[] p, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qbinom_dll(p[i], size, prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rbinom(int n, double size, double prob)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rbinom_dll(size, prob);
            return result;
        }
        #endregion Binomial Distribution----------------------------------------------------

        #region Multnomial Distribution --------------------------------------------------
        public static int[][] rmultinom(int n, int size, double[] prob)
        {
            int[][] result = new int[n][];
            for (var i = 0; i < n; i++)
            {
                //set_seed_dll((uint)(DateTime.Now.Millisecond + i), (uint)(DateTime.Now.Millisecond - i));
                int[] temp = new int[size];
                rmultinom_dll(size, prob, size, temp);
                result[i] = temp;
            }
            return result;
        }
        #endregion Multnomial Distribution --------------------------------------------------

        #region Cauchy Distribution--------------------------------------------------
        public static double dcauchy(double x, double location = 0, double scale = 1, Boolean log = false)
        {
            return dcauchy_dll(x, location, scale, log ? 1 : 0);
        }
        public static double[] dcauchy(double[] x, double location = 0, double scale = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dcauchy_dll(x[i], location, scale, log ? 1 : 0);

            return result;
        }

        public static double pcauchy(double q, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return pcauchy_dll(q, location, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pcauchy(double[] q, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pcauchy_dll(q[i], location, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qcauchy(double p, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qcauchy_dll(p, location, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qcauchy(double[] p, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qcauchy_dll(p[i], location, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rcauchy(int n, double location = 0, double scale = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rcauchy_dll(location, scale);
            return result;
        }
        #endregion Cauchy Distribution--------------------------------------------------

        #region Exponential Distribution------------------------------------------
        public static double dexp(double x, double rate = 1, Boolean log = false)
        {
            return dexp_dll(x, 1 / rate, log ? 1 : 0);
        }
        public static double[] dexp(double[] x, double rate = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dexp_dll(x[i], 1 / rate, log ? 1 : 0);

            return result;
        }

        public static double pexp(double q, double rate = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return pexp_dll(q, 1 / rate, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pexp(double[] q, double rate, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pexp_dll(q[i], 1 / rate, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qexp(double p, double rate = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qexp_dll(p, 1 / rate, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qexp(double[] p, double rate = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qexp_dll(p[i], 1 / rate, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rexp(int n, double rate = 1)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);

            for (var i = 0; i < n; i++)
                result[i] = rexp_dll(1 / rate);

            return result;
        }
        #endregion Exponential Distribution------------------------------------------

        #region Geometric Distribution---------------------------------------
        public static double dgeom(double x, double prob, Boolean log = false)
        {
            return dgeom_dll(x, prob, log ? 1 : 0);
        }
        public static double[] dgeom(double[] x, double prob, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dgeom_dll(x[i], prob, log ? 1 : 0);

            return result;
        }

        public static double pgeom(double q, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return pgeom_dll(q, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pgeom(double[] q, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pgeom_dll(q[i], prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qgeom(double p, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return qgeom_dll(p, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qgeom(double[] p, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qgeom_dll(p[i], prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rgeom(int n, double prob)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);

            for (var i = 0; i < n; i++)
                result[i] = rgeom_dll(prob);

            return result;
        }
        #endregion Geometric Distribution-------------------------------------

        #region Hypergeometric Distibution---------------------------------------------------------
        public static double dhyper(double x, double m, double n, double k, Boolean log = false)
        {
            return dhyper_dll(x, m, n, k, log ? 1 : 0);
        }
        public static double[] dhyper(double[] x, double m, double n, double k, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dhyper_dll(x[i], m, n, k, log ? 1 : 0);

            return result;
        }

        public static double phyper(double q, double m, double n, double k, Boolean lower_tail = true, Boolean log = false)
        {
            return phyper_dll(q, m, n, k, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] phyper(double[] q, double m, double n, double k, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = phyper_dll(q[i], m, n, k, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qhyper(double p, double m, double n, double k, Boolean lower_tail = true, Boolean log = false)
        {
            return qhyper_dll(p, m, n, k, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qhyper(double[] p, double m, double n, double k, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qhyper_dll(p[i], m, n, k, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rhyper(int nn, double m, double n, double k)
        {
            double[] result = new double[nn];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < nn; i++)
                result[i] = rhyper_dll(m, n, k);

            return result;
        }

        #endregion Hypergeometric Distibution---------------------------------------------------------

        #region Negative Binomial Distribution-------------------------------------------
        public static double dnbinom(double x, double size, double prob, Boolean log = false)
        {
            return dnbinom_dll(x, size, prob, log ? 1 : 0);
        }
        public static double dnbinom_mu(double x, double size, double mu, Boolean log = false)
        {
            return dnbinom_mu_dll(x, size, mu, log ? 1 : 0);
        }
        public static double[] dnbinom(double[] x, double size, double prob, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnbinom_dll(x[i], size, prob, log ? 1 : 0);

            return result;
        }
        public static double[] dnbinom_mu(double[] x, double size, double mu, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dnbinom_mu_dll(x[i], size, mu, log ? 1 : 0);

            return result;
        }

        public static double pnbinom(double q, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return pnbinom_dll(q, size, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double pnbinom_mu(double q, double size, double mu, Boolean lower_tail = true, Boolean log = false)
        {
            return pnbinom_mu_dll(q, size, mu, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pnbinom(double[] q, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnbinom_dll(q[i], size, prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] pnbinom_mu(double[] q, double size, double mu, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pnbinom_mu_dll(q[i], size, mu, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qnbinom(double p, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            return qnbinom_dll(p, size, prob, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double qnbinom_mu(double p, double size, double mu, Boolean lower_tail = true, Boolean log = false)
        {
            return qnbinom_mu_dll(p, size, mu, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qnbinom(double[] p, double size, double prob, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnbinom_dll(p[i], size, prob, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }
        public static double[] qnbinom_mu(double[] p, double size, double mu, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qnbinom_mu_dll(p[i], size, mu, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rnbinom(int n, double size, double prob)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rnbinom_dll(size, prob);

            return result;
        }

        public static double[] rnbinom_mu(int n, double size, double mu)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
            {
                result[i] = rnbinom_mu_dll(size, mu);
            }

            return result;
        }

        #endregion Negative Binomial Distribution-------------------------------------------

        #region Poisson Distribution--------------------------------------------
        public static double dpois(double x, double lambda, Boolean log = false)
        {
            return dpois_dll(x, lambda, log ? 1 : 0);
        }
        public static double[] dpois(double[] x, double lambda, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dpois_dll(x[i], lambda, log ? 1 : 0);

            return result;
        }

        public static double ppois(double q, double lambda, Boolean lower_tail = true, Boolean log = false)
        {
            return ppois_dll(q, lambda, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] ppois(double[] q, double lambda, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = ppois_dll(q[i], lambda, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qpois(double p, double lambda, Boolean lower_tail = true, Boolean log = false)
        {
            return qpois_dll(p, lambda, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qpois(double[] p, double lambda, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qpois_dll(p[i], lambda, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rpois(int n, double lambda)
        {
            double[] result = new double[n];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);

            for (var i = 0; i < n; i++)
                result[i] = rpois_dll(lambda);

            return result;
        }
        #endregion Poisson Distribution--------------------------------------------

        #region Weibull Distribution---------------------------------------
        public static double dweibull(double x, double shape, double scale = 1, Boolean log = false)
        {
            return dweibull_dll(x, shape, scale, log ? 1 : 0);
        }
        public static double[] dweibull(double[] x, double shape, double scale = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dweibull_dll(x[i], shape, scale, log ? 1 : 0);

            return result;
        }

        public static double pweibull(double q, double shape, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return pweibull_dll(q, shape, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pweibull(double[] q, double shape, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pweibull_dll(q[i], shape, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qweibull(double p, double shape, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qweibull_dll(p, shape, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qweibull(double[] p, double shape, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qweibull_dll(p[i], shape, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rweibull(int n, double shape, double scale = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rweibull_dll(shape, scale);
            return result;
        }
        #endregion Weibull Distribution---------------------------------------

        #region Logistic Distribution----------------------------------------------
        public static double dlogis(double x, double location = 0, double scale = 1, Boolean log = false)
        {
            return dlogis_dll(x, location, scale, log ? 1 : 0);
        }
        public static double[] dlogis(double[] x, double location = 0, double scale = 1, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dlogis_dll(x[i], location, scale, log ? 1 : 0);

            return result;
        }

        public static double plogis(double q, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return plogis_dll(q, location, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] plogis(double[] q, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = plogis_dll(q[i], location, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qlogis(double p, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            return qlogis_dll(p, location, scale, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qlogis(double[] p, double location = 0, double scale = 1, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qlogis_dll(p[i], location, scale, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rlogis(int n, double location = 0, double scale = 1)
        {
            double[] result = new double[n];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < n; i++)
                result[i] = rlogis_dll(location, scale);
            return result;
        }
        #endregion Logistic Distribution--------------------------------------------------

        #region Studentized Range Distribution---------------------------------
        public static double ptukey(double q, double nmeans, double df, double nranges = 1, Boolean lower_tail = true, Boolean log_p = false)
        {
            return ptukey_dll(q, nranges, nmeans, df, lower_tail ? 1 : 0, log_p ? 1 : 0);
        }
        public static double[] ptukey(double[] q, double nmeans, double df, double nranges = 1, Boolean lower_tail = true, Boolean log_p = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = ptukey_dll(q[i], nranges, nmeans, df, lower_tail ? 1 : 0, log_p ? 1 : 0);

            return result;
        }

        public static double qtukey(double p, double nmeans, double df, double nranges = 1, Boolean lower_tail = true, Boolean log_p = false)
        {
            return qtukey_dll(p, nranges, nmeans, df, lower_tail ? 1 : 0, log_p ? 1 : 0);
        }
        public static double[] qtukey(double[] p, double nmeans, double df, double nranges = 1, Boolean lower_tail = true, Boolean log_p = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qtukey_dll(p[i], nranges, nmeans, df, lower_tail ? 1 : 0, log_p ? 1 : 0);

            return result;
        }
        #endregion Studentized Range Distribution-------------------------------

        #region Wilcoxon Rank Sum Distribution-------------------------------
        public static double dwilcox(double x, double m, double n, Boolean log = false)
        {
            return dwilcox_dll(x, m, n, log ? 1 : 0);
        }
        public static double[] dwilcox(double[] x, double m, double n, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dwilcox_dll(x[i], m, n, log ? 1 : 0);

            return result;
        }

        public static double pwilcox(double q, double m, double n, Boolean lower_tail = true, Boolean log = false)
        {
            return pwilcox_dll(q, m, n, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] pwilcox(double[] q, double m, double n, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = pwilcox_dll(q[i], m, n, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qwilcox(double p, double m, double n, Boolean lower_tail = true, Boolean log = false)
        {
            return qwilcox_dll(p, m, n, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qwilcox(double[] p, double m, double n, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qwilcox_dll(p[i], m, n, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rwilcox(int nn, double m, double n)
        {
            double[] result = new double[nn];
            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);
            for (var i = 0; i < nn; i++)
                result[i] = rwilcox_dll(m, n);
            return result;
        }
        #endregion Wilcoxon Rank Sum Distribution-------------------------------

        #region Wilcoxon Signed Rank Distribution-----------------------
        public static double dsignrank(double x, double n, Boolean log = false)
        {
            return dsignrank_dll(x, n, log ? 1 : 0);
        }
        public static double[] dsignrank(double[] x, double n, Boolean log = false)
        {
            var length = x.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = dsignrank_dll(x[i], n, log ? 1 : 0);

            return result;
        }

        public static double psignrank(double q, double n, Boolean lower_tail = true, Boolean log = false)
        {
            return psignrank_dll(q, n, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] psignrank(double[] q, double n, Boolean lower_tail = true, Boolean log = false)
        {
            var length = q.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = psignrank_dll(q[i], n, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double qsignrank(double p, double n, Boolean lower_tail = true, Boolean log = false)
        {
            return qsignrank_dll(p, n, lower_tail ? 1 : 0, log ? 1 : 0);
        }
        public static double[] qsignrank(double[] p, double n, Boolean lower_tail = true, Boolean log = false)
        {
            var length = p.Length;
            double[] result = new double[length];

            for (var i = 0; i < length; i++)
                result[i] = qsignrank_dll(p[i], n, lower_tail ? 1 : 0, log ? 1 : 0);

            return result;
        }

        public static double[] rsignrank(int nn, double n)
        {
            double[] result = new double[nn];

            //set_seed_dll((uint)DateTime.Now.Millisecond, (uint)DateTime.Now.Second);

            for (var i = 0; i < nn; i++)
                result[i] = rsignrank_dll(n);

            return result;
        }
        #endregion Wilcoxon Signed Rank Distribution-----------------------

        #region Gamma and Related Functions--------------------------------------------------------------------
        public static double gamma(double x)
        {
            return gammafn_dll(x);
        }
        public static double lgamma(double x)
        {
            return lgammafn_dll(x);
        }
        public static double lgamma_sign(double x, int[] sgn)
        {
            return lgammafn_sign_dll(x, sgn);
        }

        public static double psigamma(double x, double deriv = 0)
        {
            return psigamma_dll(x, deriv);
        }
        public static double digamma(double x)
        {
            return digamma_dll(x);
        }

        public static double trigamma(double x)
        {
            return trigamma_dll(x);
        }

        public static double tetragamma(double x)
        {
            return tetragamma_dll(x);
        }

        public static double pentagamma(double x)
        {
            return pentagamma_dll(x);
        }

        public static double beta(double a, double b)
        {
            return beta_dll(a, b);
        }

        public static double lbeta(double a, double b)
        {
            return lbeta_dll(a, b);
        }

        public static double choose(double n, double k)
        {
            return choose_dll(n, k);
        }

        public static double lchoose(double n, double k)
        {
            return lchoose_dll(n, k);
        }

        #endregion Gamma and Related Functions-----------------------------------------------------------------

        #region bessel----------------------------------------------------------------------------------
        public static double besselI(double x, double nu, Boolean expon_scaled = false)
        {
            return bessel_i_dll(x, nu, 1.0 + ((expon_scaled) ? 1 : 0));
        }

        public static double besselJ(double x, double nu)
        {
            if (x < 0)
            {
                Console.Error.WriteLine("value out of range in 'bessel_j'");
                System.Environment.Exit(0);
            }
            return bessel_j_dll(x, nu);
        }

        public static double besselK(double x, double nu, Boolean expon_scaled = false)
        {
            return bessel_k_dll(x, nu, 1.0 + ((expon_scaled) ? 1 : 0));
        }

        public static double besselY(double x, double nu)
        {
            if (x < 0)
            {
                Console.Error.WriteLine("value out of range in 'bessel_y'");
                System.Environment.Exit(0);
            }
            return bessel_y_dll(x, nu);
        }

        #endregion bessel----------------------------------------------------------------------------------

        #region General Support Functions-------------------------------------------------------
        public static double hypot(double a, double b)
        {
            return hypot_dll(a, b);
        }

        public static double pythag(double a, double b)
        {
            return pythag_dll(a, b);
        }

        public static double expm1(double x)
        {
            return expm1_dll(x);
        }

        public static double log1p(double x)
        {
            return log1p_dll(x);
        }

        public static int imax2(int a, int b)
        {
            return imax2_dll(a, b);
        }

        public static int imin2(int a, int b)
        {
            return imin2_dll(a, b);
        }

        public static double fmax2(double a, double b)
        {
            return fmax2_dll(a, b);
        }

        public static double fmin2(double a, double b)
        {
            return fmin2_dll(a, b);
        }

        public static double sign(double x)
        {
            return sign_dll(x);
        }

        public static double fprec(double a, double b)
        {
            return fprec_dll(a, b);
        }

        public static double round(double a, double b)
        {
            return fround_dll(a, b);
        }

        public static double fsign(double a, double b)
        {
            return fsign_dll(a, b);
        }

        public static double trunc(double x)
        {
            return ftrunc_dll(x);
        }

        public static double log1pmx(double x)
        {
            return log1pmx_dll(x);
        }

        public static double lgamma1p(double x)
        {
            return lgamma1p_dll(x);
        }

        public static double logspace_add(double a, double b)
        {
            return logspace_add_dll(a, b);
        }

        public static double logspace_sub(double a, double b)
        {
            return logspace_sub_dll(a, b);
        }

        #endregion General Support Functions------------------------------------------------------------------

        #endregion wrapper------------------------------------------------------------------------------------
    }
}
